1. Open https://github.com/dlab-geo/r-geospatial-workshop
   • Download & unzip the zip file
   • Make a note of the folder in which it is located

2. Start RStudio and open a new script file

3. Install required libraries in RStudio - if you do not have them already

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

Intro to working with geospatial data in R

• geospatial data files and formats
• Loading geospatial data in R
• R packages for working with geospatial data
• coordinate reference systems

Mapping geospatial data

Practice

Patty Frontiera

dlab.berkeley.edu

Who are you?

Why are you here?

Open one of the tutorial files in a web browser

Slides r-geospatial-workshop-pt1.html

Tutorial r-geospatial-workshop-pt1-tutorial.html

Raw Code scripts/r-geospatial-workshop-pt1.R

Make sure you can cut and paste into RStudio # Geographic Data

are data about locations on or near the surface of the Earth.

represent location more specifically with coordinates

46.130479, -117.134167

Coordinates only make sense when associated with a CRS!

Geographic Coordinates: Latitude and Longitude

Define:

• the shape of the Earth

• the origin (0,0 point)

• the relationship between the system and the real world

• the units

Because of variations in these, there are many geographic CRSs!

The World Geodetic System of 1984 is the most widely used geographic coordinate reference system.

WGS84 is the default CRS for most GIS software

Almost all longitude and latitude data are assumed to be WGS84 unless otherwise specified

Historical data much trickier

You can

• dynamically determine spatial metrics like area, length, distance and direction

• spatial relationships like intersects, inside, contains, etc

and

• link data by location, like census data and crime data

Spatial data is a broader term than geographic data.

Methods for working with spatial data are valid for geospatial data

All spatial data are encoded with some type of coordinate reference system

Geospatial data require the CRS to be a model of locations on the Earth

Vector and Raster Data

Points, lines and Polygons

Regular grid of cells (or pixels)

• We won't be covering Raster data in this workshop*

software that can import, create, store, edit, visualize and analyze geospatial data

• represented as geometric data objects referenced to the surface of the earth via CRSs

• with methods to operate on those representations

We call software for working with geospatial data GIS

Geographic Information System

This term is commonly associated with desktop software applications.

Desktop GIS - ArcGIS, QGIS

Spatial Databases - PostgreSQL/PostGIS

Web-based GIS - ArcGIS Online, CARTO

Software geospatial data support - Tableau

Programming languages with geospatial data support

• R, Python, Javascript

You already use R

Reproducibility

Free & Open Source

Strong support for geospatial data and analysis

Cutting edge

Make sure you have the packages we are going to use installed and loaded

Make sure you have set your working directory to the location of the workshop files

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

There are many approaches to and packages for working with geospatial data in R.

One approach is to keep it simple and store geospatial data in a data frame.

This approach is most common when

• the data are point data in CSV files and

• you want to map rather than spatially transform or analyze the data

sf_properties_25ksample.csv

San Francisco Open Data Portal https://data.sfgov.org

SF Property Tax Rolls

This data set includes the Office of the Assessor-Recorder’s secured property tax roll spanning from 2007 to 2016.

We are using a subset of these data as a proxy for home values.

sfhomes <- read.csv('data/sf_properties_25ksample.csv', 
                    stringsAsFactors = FALSE)

# Take a look at first 5 rows and a few of the columns
sfhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
              "NumBedrooms")]

Make sure your working directory is set to the folder where you downloaded the workshop files!

sfhomes <- read.csv('data/sf_properties_25ksample.csv', 
                    stringsAsFactors = FALSE)

# Take a look at first 5 rows and a few of the columns
sfhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
              "NumBedrooms")]

##   YearBuilt totvalue AreaSquareFeet       Neighborhood NumBedrooms
## 1      1923  1031391           2250 West of Twin Peaks           0
## 2      1909   775300           5830   Presidio Heights           3
## 3      1915  1116772           1350     Inner Richmond           2
## 4      1947   860130           1162    Sunset/Parkside           0
## 5      1981   322749           1463     Outer Richmond           3

class(sfhomes)            # what is the data object type?
dim(sfhomes)              # how many rows and columns
str(sfhomes)              # display the structure of the object
head(sfhomes)             # take a look at the first 10 records
summary(sfhomes)          # explore the range of values
summary(sfhomes$totvalue) # explore the range of values for one column
hist(sfhomes$totvalue)    # histogram for the totvalue column

• What columns contain the geographic data?
• Are these data vector or raster data?
• What type of geometry do the data contain?
  • Points, lines, polygons, grid cells?
• What is the CRS of these data?

Use the R base plot function to create a simple map

plot(sfhomes$lon, sfhomes$lat) # using base plot function

Use the R base plot function to create a simple map

plot(sfhomes$lon, sfhomes$lat) # using base plot function

Most widely used plotting library in R

Not specifically for geospatial data

But can be used to make fabulous maps

Great choice if you already know ggplot2

Load the library

library(ggplot2)

Basic map with ggplot

library(ggplot2)

ggplot() + geom_point(data=sfhomes, aes(lon,lat))

Basic map with ggplot

ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1)

ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1) + coord_map()

Allows you to associate a map projection with geographic coord data.

Map Projection: mathematial transformation from curved to flat surface

A Projected CRS applies a map projection to a Geographic CRS

All introduce distortion,

• in shape, area, distance, direction, or combo

• the larger the area the greater the distortion

No one map projection best for all purposes

Selection depends on location, extent and purpose

Data driven symbology

ggplot() + geom_point(data=sfhomes, aes(lon,lat, col=totvalue)) + 
  coord_map()

What's happening here?

sfhomes_low2high <- sfhomes[order(sfhomes$totvalue, decreasing = FALSE),]

ggplot() + 
  geom_point(data=sfhomes_low2high, aes(lon,lat, col=totvalue)) + 
  coord_map()

Try it - Does the output map look different from previous one?

The order of the data in the data frame changes the map display!

Map the sfhomes data in decreasing order by totvalue.

sfhomes_high2low <- sfhomes[order(sfhomes$totvalue, decreasing = T),]
ggplot() + geom_point(data=sfhomes_high2low, aes(lon,lat, col=totvalue)) + 
  coord_map()

What does this code do?

sfhomes2010_15 <- subset(sfhomes_low2high, as.numeric(SalesYear) > 2009)

ggplot() +
  geom_point(aes(lon, lat, col=totvalue), data = sfhomes2010_15 )  +
  facet_wrap(~ SalesYear)

Visual spatial analysis!

Create basemaps on which you can display your data.

Geocode place names and addresses to get point coordinates.

and more…

Load the libary

library(ggmap)

## Google Maps API Terms of Service: http://developers.google.com/maps/terms.

## Please cite ggmap if you use it: see citation("ggmap") for details.

Some ggmap functionality may require you to register a Google API key

• Search - Google API Key
• Create Google Maps API Key
• Enable API key for geocoding, static maps API, javascript API

with ggmap

#register_google(key="AIzXXXXXXXXXXXXXXXXXxPSE") # your key here
register_google(key="AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo")

This key will be disabled after the workshop!

If you have problems with ggmap you may need to reinstall the library from github.

See this StackOverflow discussion

#devtools::install_github("dkahle/ggmap")
#library(ggmap)

Fetch map data to plot using get_map

Use ?get_map for details

sf_map <- get_map("San Francisco, CA")  

This use of the command leverages the Google Geocoding API to determine the coordinates of the named location.

Fetch map data to plot using get_map

Use ?get_map for details

sf_map <- get_map("San Francisco, CA")  

## Source : https://maps.googleapis.com/maps/api/staticmap?center=San+Francisco,+CA&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=San%20Francisco%2C%20CA

ggmap(sf_map)

ggmap(sf_map)

Syntax similar to ggplot

# ggplot() +

ggmap(sf_map) + 
  geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))

ggmap(sf_map) +
  geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))

Create a basemap zoomed to the extent of our data

# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes$lon), lat = mean(sfhomes$lat))
sf_ctr  # take a look

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

Get a basemap zoomed to the extent of the 2015 data

# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes15$lon), lat = mean(sfhomes15$lat))
sf_ctr  # take a look

##        lon        lat 
## -122.43289   37.75964

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=37.759644,-122.432893&zoom=12&size=640x640&scale=1&maptype=terrain&language=en-EN

Create a ggmap that uses the customized sf_basemap

• Recreate our last map of property values (totvalue)

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

Finally, let's add another geospatial data layer to our ggmap.

You can use this method to add as many layers as you want to a ggmap or ggplot.

Use the read.csv function to read in a file of Bart Station locations.

bart <- read.csv("./data/bart.csv")
# take a look
head (bart)

##           X        Y             STATION OPERATOR DIST  CO
## 1 -122.2833 37.87406      NORTH BERKELEY     BART    4 ALA
## 2 -122.2682 37.86969   DOWNTOWN BERKELEY     BART    4 ALA
## 3 -122.2701 37.85321               ASHBY     BART    4 ALA
## 4 -122.2518 37.84451           ROCKRIDGE     BART    4 ALA
## 5 -122.2671 37.82871           MACARTHUR     BART    4 ALA
## 6 -122.2684 37.80808 19TH STREET/OAKLAND     BART    4 ALA

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red") 

## Warning: Removed 34 rows containing missing values (geom_point).

Questons?

Determining the geographic coordinates for one or more addresses, zip codes, or named places.

police_stations <- read.csv("data/sf_police_addresses.csv", 
                            stringsAsFactors = F)
head(police_stations)

##   PoliceDistrict        Phone              Address          City State
## 1        Central 415-315-2400   766 Vallejo Street San Francisco    CA
## 2       Southern 415-575-6000      1251 3rd Street San Francisco    CA
## 3        Bayview 415-671-2300  201 Williams Avenue San Francisco    CA
## 4        Mission 415-558-5400  630 Valencia Street San Francisco    CA
## 5       Northern 415-614-3400 1125 Fillmore Street San Francisco    CA
## 6           Park 415-242-3000   1899 Waller Street San Francisco    CA

#?geocode
geocode(police_stations)

police_stations$full_address <- paste(police_stations$Address, 
                                      police_stations$City, police_stations$State)
head(police_stations)

##   PoliceDistrict        Phone              Address          City State
## 1        Central 415-315-2400   766 Vallejo Street San Francisco    CA
## 2       Southern 415-575-6000      1251 3rd Street San Francisco    CA
## 3        Bayview 415-671-2300  201 Williams Avenue San Francisco    CA
## 4        Mission 415-558-5400  630 Valencia Street San Francisco    CA
## 5       Northern 415-614-3400 1125 Fillmore Street San Francisco    CA
## 6           Park 415-242-3000   1899 Waller Street San Francisco    CA
##                            full_address
## 1   766 Vallejo Street San Francisco CA
## 2      1251 3rd Street San Francisco CA
## 3  201 Williams Avenue San Francisco CA
## 4  630 Valencia Street San Francisco CA
## 5 1125 Fillmore Street San Francisco CA
## 6   1899 Waller Street San Francisco CA

station_coords <- geocode(police_stations$full_address)

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=766%20Vallejo%20Street%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1251%203rd%20Street%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=201%20Williams%20Avenue%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=630%20Valencia%20Street%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1125%20Fillmore%20Street%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1899%20Waller%20Street%20San%20Francisco%20CA

## Warning: geocode failed with status OVER_QUERY_LIMIT, location = "1899
## Waller Street San Francisco CA"

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=461%206th%20Avenue%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1%20Sgt.%20John%20V.%20Young%20Lane%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=2345%2024th%20Avenue%20San%20Francisco%20CA

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=301%20Eddy%20Street%20San%20Francisco%20CA

## Warning: geocode failed with status OVER_QUERY_LIMIT, location = "301 Eddy
## Street San Francisco CA"

police_stations <- cbind(police_stations, station_coords)
police_stations[1:5,c("PoliceDistrict","full_address","lon","lat")]

##   PoliceDistrict                          full_address       lon      lat
## 1        Central   766 Vallejo Street San Francisco CA -122.4100 37.79872
## 2       Southern      1251 3rd Street San Francisco CA -122.3894 37.77238
## 3        Bayview  201 Williams Avenue San Francisco CA -122.3980 37.72976
## 4        Mission  630 Valencia Street San Francisco CA -122.4220 37.76285
## 5       Northern 1125 Fillmore Street San Francisco CA -122.4325 37.78022

to your last map

Use a different color

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red", size=3) +
  geom_point(data=police_stations, aes(x=lon,y=lat), shape=22, 
             col="black", fill="grey", size=4)

## Warning: Removed 34 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

## About Google geocoding

• You can geocode addresses or place names

• Super accurate

• You may need an API key

• Limited to 2500 per day

• Use ESRI Geocoding tools if you have more addresses

SF Landmarks

data/landmarks.csv

landmarks <- read.csv("./data/landmarks.csv")
head(landmarks)

##           X       Y id             name
## 1 -13626151 4551548  1       Coit Tower
## 2 -13630804 4544523  2       Twin Peaks
## 3 -13627654 4548289  3        City Hall
## 4 -13635101 4546904  4 Golden Gate Park

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red", size=3) +
  geom_point(data=landmarks, aes(x=X,y=Y), shape=22, 
             col="black", fill="grey", size=4)

All good? If not, why not?

## Warning: Removed 34 rows containing missing values (geom_point).

## Warning: Removed 4 rows containing missing values (geom_point).

There are limits to what you can do with geospatial data stored in a data frame.

These data are not in geographic coordinates - longitude and latitude.

You can't map these with ggmap.

##           X       Y id             name
## 1 -13626151 4551548  1       Coit Tower
## 2 -13630804 4544523  2       Twin Peaks
## 3 -13627654 4548289  3        City Hall
## 4 -13635101 4546904  4 Golden Gate Park

The ESRI Shapefile is the most common file format for geospatial data.

Can't read in or directly plot data in a shapefile with ggplot or ggmap.

Spatial data objects and methods are needed to answer questions like

• What properties are in the Noe Valley neighborhood?

• What is the average property value in each SF neighborhood?

• What is the area of each SF Neighborhood and the property density?

• What properties are within walking distance (.25 miles) of the Mission neighborhood?

Classes and Methods for Spatial Data

The SP package is most commonly used to construct and manipulate spatial data objects in R.

Hundreds of other R packages that do things with spatial data typically build on SP objects.

Load the library sp

Take a look at the different types of spatial object classes supported by sp

library(sp)
getClass("Spatial") 

## Warning: package 'sp' was built under R version 3.4.3

## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

We use the S*DF objects most frequently!

Let's transform the sfhomes15 data frame to an sp object of type SpatialPointsDataFrame

Use the sp::coordinates() method

• Sets or retrieves spatial coordinates

When transforming a DF to SPDF you are setting the coordinates.

Note - the transformation happens in place!

sfhomes15_sp <- sfhomes15  # Make a copy - why?

class(sfhomes15_sp)  # check the class of the object

coordinates(sfhomes15_sp) <- c('lon','lat')  #  Make it spatial - ORDER MATTERS!!

class(sfhomes15_sp)  # check the class of the object

sfhomes15_sp <- sfhomes15  # Make a copy - why?

class(sfhomes15_sp)  # check the class of the object

## [1] "data.frame"

coordinates(sfhomes15_sp) <- c('lon','lat')  #  Make it spatial - ORDER MATTERS!!

class(sfhomes15_sp)  # check the class of the object

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

You transformed a data frame to an SPDF using:

coordinates(sfhomes15_sp) <- c('lon','lat')

Now try this:

coordinates(sfhomes15_sp)

coordinates(sfhomes15_sp)

##             lon      lat
## 24    -122.4741 37.73601
## 35    -122.3965 37.75920
## 76    -122.3942 37.77181
## 79    -122.4754 37.71518
## 85    -122.4210 37.75876
## 90    -122.4149 37.79362
## 92    -122.4370 37.74020
## 98    -122.4437 37.78881
## 223   -122.4754 37.71492
## 229   -122.3919 37.78868
## 273   -122.3968 37.77945
## 341   -122.5040 37.75515
## 363   -122.4041 37.76422
## 366   -122.3934 37.77764
## 388   -122.4166 37.74931
## 451   -122.4112 37.73683
## 467   -122.4507 37.76448
## 473   -122.3990 37.75450
## 487   -122.4195 37.76649
## 534   -122.4225 37.78676
## 548   -122.4410 37.70854
## 549   -122.5009 37.74485
## 550   -122.4220 37.73048
## 559   -122.4947 37.75832
## 567   -122.3942 37.77181
## 595   -122.4187 37.80469
## 624   -122.4019 37.77323
## 631   -122.4229 37.80586
## 646   -122.4225 37.78676
## 658   -122.3919 37.78868
## 665   -122.4261 37.74349
## 668   -122.3904 37.78139
## 683   -122.4458 37.74398
## 774   -122.3883 37.73754
## 776   -122.4767 37.75347
## 864   -122.4326 37.74229
## 974   -122.4380 37.72624
## 1018  -122.4452 37.76906
## 1097  -122.4937 37.77999
## 1155  -122.5047 37.75390
## 1164  -122.4498 37.75451
## 1167  -122.4055 37.80424
## 1188  -122.4194 37.75596
## 1242  -122.4088 37.77953
## 1273  -122.4331 37.76019
## 1374  -122.3918 37.77899
## 1379  -122.4070 37.75674
## 1487  -122.4236 37.77800
## 1598  -122.4758 37.75497
## 1600  -122.4419 37.74582
## 1618  -122.5030 37.74587
## 1668  -122.3989 37.78990
## 1698  -122.4236 37.77800
## 1747  -122.3972 37.76314
## 1755  -122.4434 37.75159
## 1768  -122.4242 37.80551
## 1804  -122.4549 37.76402
## 1828  -122.4320 37.78195
## 1833  -122.4405 37.71825
## 1913  -122.4146 37.74714
## 1999  -122.3942 37.77181
## 2024  -122.4676 37.71048
## 2033  -122.4249 37.76548
## 2066  -122.4243 37.75206
## 2073  -122.4097 37.77613
## 2096  -122.3996 37.73859
## 2248  -122.4045 37.79929
## 2258  -122.4442 37.78322
## 2285  -122.3919 37.78868
## 2314  -122.4026 37.77879
## 2329  -122.4368 37.78228
## 2339  -122.4240 37.74366
## 2355  -122.4126 37.79035
## 2416  -122.4579 37.72632
## 2428  -122.4089 37.79286
## 2436  -122.4409 37.80451
## 2443  -122.4127 37.78995
## 2500  -122.4357 37.76468
## 2508  -122.4666 37.76347
## 2537  -122.4870 37.74605
## 2538  -122.3942 37.77181
## 2627  -122.4219 37.77389
## 2692  -122.3958 37.77705
## 2717  -122.4116 37.77806
## 2734  -122.4290 37.72977
## 2793  -122.4764 37.73390
## 2914  -122.4183 37.76712
## 2969  -122.4711 37.72241
## 2975  -122.4393 37.78018
## 3127  -122.4393 37.78951
## 3157  -122.3702 37.72838
## 3205  -122.4398 37.78746
## 3262  -122.4540 37.78251
## 3288  -122.4698 37.75298
## 3296  -122.4607 37.73104
## 3306  -122.4081 37.77962
## 3354  -122.4849 37.74663
## 3403  -122.4209 37.75792
## 3436  -122.4854 37.74788
## 3471  -122.4856 37.78116
## 3485  -122.4177 37.80345
## 3516  -122.4354 37.72201
## 3578  -122.4199 37.77606
## 3611  -122.4269 37.80229
## 3658  -122.4108 37.76534
## 3754  -122.3888 37.72872
## 3772  -122.3891 37.73264
## 4023  -122.4827 37.77762
## 4025  -122.4450 37.76430
## 4058  -122.4328 37.71546
## 4076  -122.4660 37.74828
## 4111  -122.5023 37.74333
## 4114  -122.4970 37.74468
## 4140  -122.4741 37.73601
## 4150  -122.4060 37.79853
## 4249  -122.4302 37.74726
## 4270  -122.4194 37.76710
## 4287  -122.4676 37.71048
## 4369  -122.4950 37.73327
## 4376  -122.3919 37.78868
## 4393  -122.3972 37.77570
## 4427  -122.4332 37.77869
## 4466  -122.4786 37.77394
## 4469  -122.4210 37.76470
## 4483  -122.4318 37.71761
## 4514  -122.4214 37.78139
## 4526  -122.4153 37.74047
## 4527  -122.4209 37.79368
## 4555  -122.4145 37.79167
## 4561  -122.3910 37.78938
## 4566  -122.4031 37.79058
## 4578  -122.3906 37.75501
## 4579  -122.4154 37.74537
## 4592  -122.3922 37.78607
## 4599  -122.4385 37.75980
## 4602  -122.4785 37.76295
## 4618  -122.4137 37.79099
## 4648  -122.4635 37.76461
## 4663  -122.4160 37.74563
## 4669  -122.3918 37.77899
## 4759  -122.4911 37.75318
## 4800  -122.4137 37.79099
## 4825  -122.4147 37.72528
## 4853  -122.4388 37.75978
## 4864  -122.3696 37.72900
## 4896  -122.4245 37.79493
## 4899  -122.3902 37.75236
## 4921  -122.5103 37.77369
## 4939  -122.3930 37.73299
## 4944  -122.3919 37.78868
## 4970  -122.4646 37.77383
## 4986  -122.3904 37.78139
## 5053  -122.4287 37.77484
## 5085  -122.4244 37.78587
## 5104  -122.4255 37.72744
## 5108  -122.4364 37.77781
## 5109  -122.4394 37.80326
## 5114  -122.4853 37.75287
## 5188  -122.4302 37.74806
## 5228  -122.3881 37.73174
## 5247  -122.4242 37.80551
## 5252  -122.4296 37.79988
## 5261  -122.4048 37.79879
## 5287  -122.4335 37.79946
## 5288  -122.4611 37.73481
## 5316  -122.3918 37.77899
## 5320  -122.4481 37.77945
## 5321  -122.4477 37.71766
## 5334  -122.3919 37.78868
## 5350  -122.4336 37.72535
## 5351  -122.4540 37.78251
## 5360  -122.4264 37.72624
## 5366  -122.3922 37.78607
## 5398  -122.3905 37.78759
## 5425  -122.3919 37.78868
## 5430  -122.4263 37.72116
## 5450  -122.4175 37.79865
## 5461  -122.3989 37.78577
## 5467  -122.4175 37.77341
## 5511  -122.3919 37.78868
## 5519  -122.4093 37.74024
## 5521  -122.4816 37.77294
## 5522  -122.4194 37.75596
## 5547  -122.3959 37.79054
## 5561  -122.4213 37.78884
## 5571  -122.4194 37.76710
## 5609  -122.4245 37.79920
## 5647  -122.4085 37.75911
## 5654  -122.4673 37.75439
## 5671  -122.4459 37.75811
## 5691  -122.4806 37.77724
## 5705  -122.4167 37.76091
## 5712  -122.3702 37.72838
## 5732  -122.4207 37.73596
## 5740  -122.4101 37.77808
## 5757  -122.4752 37.71513
## 5765  -122.3922 37.73203
## 5827  -122.4036 37.71724
## 5978  -122.4264 37.73053
## 5997  -122.3886 37.77159
## 6043  -122.5068 37.77785
## 6056  -122.4089 37.79286
## 6134  -122.4191 37.70898
## 6140  -122.4107 37.76483
## 6181  -122.3974 37.78575
## 6191  -122.4321 37.73585
## 6208  -122.4774 37.75583
## 6252  -122.4355 37.74743
## 6258  -122.3904 37.78271
## 6342  -122.4846 37.75914
## 6365  -122.4356 37.75016
## 6410  -122.3824 37.72887
## 6430  -122.4699 37.75813
## 6477  -122.4656 37.71671
## 6581  -122.4792 37.77343
## 6687  -122.4243 37.79352
## 6691  -122.4371 37.71173
## 6694  -122.4399 37.71837
## 6721  -122.4308 37.73647
## 6723  -122.4566 37.75242
## 6782  -122.4298 37.71104
## 6787  -122.4405 37.80504
## 6814  -122.4848 37.73940
## 6863  -122.4122 37.79293
## 6892  -122.4419 37.74582
## 6903  -122.4084 37.79154
## 6915  -122.4144 37.77516
## 6918  -122.3958 37.77816
## 6932  -122.4244 37.74805
## 6933  -122.4368 37.78249
## 6988  -122.5064 37.77555
## 7006  -122.4528 37.76708
## 7009  -122.4947 37.76067
## 7032  -122.4050 37.77762
## 7073  -122.4194 37.75596
## 7142  -122.4953 37.73096
## 7149  -122.4104 37.75357
## 7175  -122.4370 37.76934
## 7206  -122.4936 37.77938
## 7303  -122.4327 37.77299
## 7320  -122.4194 37.76710
## 7337  -122.4752 37.71551
## 7395  -122.4263 37.76856
## 7404  -122.3939 37.75915
## 7447  -122.4113 37.72850
## 7462  -122.4127 37.76517
## 7493  -122.4194 37.75596
## 7533  -122.4515 37.75418
## 7555  -122.4071 37.77553
## 7564  -122.4389 37.73129
## 7600  -122.4097 37.77613
## 7649  -122.4290 37.78872
## 7651  -122.3958 37.77816
## 7704  -122.4337 37.76214
## 7735  -122.3904 37.78271
## 7832  -122.4853 37.75355
## 7833  -122.3983 37.75753
## 7848  -122.4490 37.76795
## 7936  -122.4087 37.77310
## 7981  -122.4489 37.71967
## 7991  -122.4325 37.77805
## 8047  -122.4019 37.78005
## 8089  -122.4194 37.75596
## 8127  -122.4685 37.76188
## 8161  -122.4307 37.74163
## 8275  -122.4126 37.77870
## 8322  -122.3924 37.73700
## 8344  -122.4126 37.77870
## 8378  -122.4900 37.78324
## 8401  -122.4642 37.71638
## 8410  -122.4279 37.77408
## 8419  -122.4016 37.75145
## 8422  -122.4836 37.71006
## 8438  -122.4668 37.73426
## 8517  -122.3919 37.78868
## 8537  -122.3980 37.71952
## 8590  -122.4291 37.73289
## 8612  -122.4194 37.75596
## 8652  -122.4222 37.72570
## 8694  -122.3919 37.78868
## 8703  -122.3958 37.77339
## 8713  -122.4301 37.77286
## 8730  -122.4126 37.77870
## 8818  -122.3897 37.78278
## 8823  -122.4361 37.78937
## 8842  -122.4257 37.79707
## 8844  -122.4259 37.78673
## 8888  -122.4318 37.79685
## 8891  -122.4277 37.77402
## 8899  -122.3896 37.78063
## 8907  -122.4728 37.75198
## 8932  -122.4102 37.74514
## 8962  -122.3942 37.77181
## 8983  -122.4206 37.80236
## 9049  -122.4651 37.75119
## 9054  -122.4371 37.73361
## 9179  -122.4262 37.72086
## 9216  -122.4160 37.73419
## 9219  -122.4550 37.74285
## 9356  -122.4522 37.73356
## 9376  -122.3982 37.75309
## 9384  -122.4152 37.73405
## 9469  -122.4174 37.80442
## 9479  -122.4251 37.72966
## 9490  -122.4835 37.74208
## 9569  -122.4723 37.73269
## 9639  -122.4883 37.77411
## 9640  -122.4509 37.73126
## 9647  -122.3909 37.75967
## 9675  -122.3891 37.76366
## 9691  -122.4381 37.78210
## 9764  -122.4875 37.76255
## 9776  -122.4422 37.72376
## 9806  -122.3962 37.77393
## 9830  -122.4642 37.71673
## 9906  -122.4433 37.76135
## 9959  -122.3934 37.77764
## 9966  -122.4639 37.72971
## 9985  -122.4411 37.75920
## 10010 -122.4279 37.76639
## 10046 -122.4755 37.71495
## 10063 -122.3879 37.73919
## 10066 -122.4429 37.78894
## 10100 -122.3890 37.73260
## 10157 -122.4377 37.79704
## 10164 -122.4180 37.74667
## 10170 -122.4453 37.79137
## 10229 -122.4137 37.79099
## 10239 -122.4451 37.78484
## 10257 -122.4102 37.74590
## 10263 -122.4353 37.75147
## 10272 -122.4974 37.78118
## 10273 -122.4401 37.75541
## 10307 -122.4167 37.79138
## 10338 -122.4343 37.73743
## 10343 -122.4108 37.77597
## 10345 -122.5037 37.74846
## 10351 -122.3910 37.78938
## 10353 -122.4865 37.76343
## 10410 -122.3940 37.75300
## 10415 -122.4287 37.73629
## 10419 -122.4344 37.76357
## 10499 -122.4531 37.73537
## 10530 -122.4150 37.73959
## 10578 -122.4594 37.73954
## 10583 -122.4760 37.71558
## 10600 -122.4832 37.78375
## 10607 -122.4911 37.73160
## 10632 -122.4347 37.78713
## 10646 -122.4518 37.76347
## 10661 -122.4066 37.73036
## 10781 -122.4124 37.80434
## 10782 -122.4554 37.71773
## 10800 -122.4152 37.73705
## 10857 -122.4178 37.76399
## 10887 -122.4260 37.74559
## 10914 -122.4349 37.74579
## 10957 -122.4173 37.74421
## 10967 -122.4587 37.73858
## 10985 -122.4325 37.77540
## 11005 -122.4558 37.73265
## 11023 -122.4173 37.73753
## 11035 -122.4386 37.76810
## 11094 -122.4275 37.72671
## 11158 -122.3937 37.78702
## 11160 -122.4236 37.77800
## 11170 -122.5054 37.76363
## 11171 -122.4001 37.79803
## 11180 -122.4757 37.71556
## 11181 -122.4359 37.75531
## 11202 -122.4133 37.78070
## 11258 -122.4437 37.80053
## 11265 -122.4865 37.76465
## 11284 -122.4117 37.75704
## 11286 -122.4194 37.75596
## 11305 -122.4136 37.79298
## 11385 -122.4424 37.74562
## 11451 -122.4370 37.73761
## 11498 -122.4438 37.71852
## 11554 -122.4990 37.77913
## 11564 -122.4194 37.75596
## 11576 -122.3910 37.78938
## 11598 -122.3824 37.72852
## 11603 -122.4748 37.71490
## 11622 -122.3988 37.71372
## 11636 -122.4678 37.72225
## 11690 -122.4415 37.78044
## 11728 -122.4793 37.75869
## 11729 -122.4165 37.73373
## 11808 -122.4457 37.75769
## 11830 -122.4965 37.77600
## 11844 -122.4468 37.73736
## 11873 -122.4474 37.73477
## 11876 -122.4126 37.79035
## 11882 -122.3970 37.75975
## 11918 -122.4438 37.76155
## 11927 -122.4333 37.73432
## 11968 -122.3936 37.78370
## 11969 -122.3907 37.78117
## 11971 -122.4742 37.75636
## 11972 -122.4542 37.77443
## 12052 -122.4385 37.75442
## 12232 -122.5060 37.73665
## 12285 -122.4786 37.77955
## 12350 -122.4672 37.73270
## 12383 -122.4093 37.77680
## 12428 -122.4274 37.77132
## 12467 -122.3971 37.78091
## 12533 -122.4024 37.77546
## 12559 -122.4130 37.80397
## 12569 -122.4114 37.72936
## 12622 -122.5059 37.73682
## 12668 -122.4131 37.79777
## 12672 -122.4154 37.73877
## 12675 -122.4357 37.73740
## 12730 -122.4011 37.73263
## 12770 -122.4221 37.72534
## 12783 -122.3910 37.78938
## 12865 -122.3919 37.78868
## 12866 -122.4891 37.75507
## 13000 -122.3968 37.77945
## 13015 -122.4341 37.76397
## 13054 -122.4165 37.77528
## 13177 -122.4043 37.78622
## 13251 -122.4294 37.77067
## 13252 -122.3933 37.77891
## 13268 -122.4316 37.76576
## 13293 -122.4936 37.77938
## 13324 -122.4116 37.75238
## 13382 -122.4513 37.76809
## 13391 -122.4248 37.80057
## 13414 -122.4433 37.78314
## 13422 -122.4141 37.80393
## 13470 -122.4439 37.71527
## 13472 -122.4815 37.74502
## 13478 -122.4214 37.78139
## 13507 -122.4298 37.73805
## 13529 -122.4242 37.80551
## 13550 -122.4439 37.80384
## 13581 -122.3956 37.73468
## 13593 -122.4286 37.74623
## 13595 -122.4617 37.78577
## 13662 -122.3918 37.77899
## 13675 -122.3696 37.72900
## 13701 -122.3962 37.77393
## 13722 -122.4025 37.73527
## 13777 -122.4750 37.71499
## 13821 -122.4032 37.78824
## 13895 -122.4485 37.73998
## 13903 -122.3919 37.78868
## 13926 -122.4141 37.75915
## 13936 -122.4296 37.73278
## 14044 -122.4484 37.76533
## 14072 -122.4546 37.74500
## 14087 -122.4642 37.75929
## 14239 -122.4665 37.71392
## 14241 -122.4456 37.80221
## 14256 -122.4990 37.77356
## 14287 -122.4234 37.77253
## 14338 -122.4475 37.77957
## 14345 -122.4192 37.75400
## 14353 -122.3696 37.72900
## 14367 -122.4253 37.79247
## 14413 -122.4397 37.76158
## 14496 -122.4276 37.79552
## 14517 -122.4617 37.75037
## 14679 -122.4565 37.72132
## 14690 -122.4365 37.71167
## 14699 -122.4375 37.76305
## 14721 -122.3919 37.78868
## 14728 -122.3934 37.77764
## 14769 -122.4258 37.78636
## 14773 -122.4326 37.75640
## 14777 -122.3936 37.78370
## 14803 -122.4121 37.77255
## 14870 -122.4481 37.73958
## 14965 -122.4093 37.74690
## 14974 -122.5070 37.78027
## 15069 -122.4340 37.74511
## 15076 -122.4733 37.77618
## 15143 -122.4234 37.77253
## 15159 -122.4937 37.78119
## 15167 -122.4129 37.78791
## 15179 -122.4668 37.74374
## 15186 -122.4476 37.77794
## 15187 -122.3991 37.78318
## 15197 -122.3943 37.77898
## 15208 -122.4557 37.71642
## 15217 -122.4678 37.72282
## 15260 -122.4416 37.71429
## 15267 -122.4507 37.74389
## 15319 -122.4345 37.78565
## 15327 -122.4957 37.77863
## 15328 -122.4564 37.70934
## 15349 -122.4151 37.73756
## 15357 -122.3761 37.73193
## 15426 -122.4009 37.76461
## 15474 -122.4009 37.76461
## 15485 -122.5046 37.75155
## 15486 -122.4718 37.76104
## 15498 -122.4174 37.74938
## 15525 -122.4124 37.77093
## 15551 -122.3943 37.78613
## 15567 -122.4055 37.80424
## 15580 -122.3914 37.78451
## 15587 -122.5046 37.74607
## 15589 -122.4339 37.72119
## 15604 -122.4248 37.75359
## 15615 -122.4098 37.71496
## 15621 -122.4667 37.73841
## 15666 -122.4379 37.75930
## 15669 -122.4194 37.75596
## 15713 -122.3702 37.72838
## 15717 -122.4382 37.75646
## 15790 -122.4021 37.76545
## 15887 -122.4311 37.76352
## 15918 -122.3942 37.77181
## 15933 -122.4178 37.76399
## 15945 -122.4390 37.79617
## 16033 -122.4014 37.78629
## 16076 -122.4149 37.79362
## 16105 -122.4489 37.76179
## 16173 -122.4135 37.72849
## 16199 -122.4014 37.77974
## 16204 -122.4327 37.76614
## 16210 -122.4210 37.75876
## 16212 -122.4126 37.77870
## 16215 -122.3958 37.77705
## 16229 -122.4400 37.71831
## 16310 -122.3988 37.72731
## 16315 -122.4254 37.80339
## 16389 -122.4800 37.75339
## 16440 -122.4970 37.74495
## 16450 -122.4194 37.75596
## 16458 -122.4624 37.77409
## 16474 -122.4542 37.77443
## 16480 -122.4014 37.77974
## 16500 -122.4126 37.77870
## 16507 -122.4300 37.78804
## 16517 -122.4240 37.76628
## 16531 -122.4266 37.73311
## 16549 -122.5038 37.74097
## 16640 -122.4384 37.77517
## 16643 -122.4517 37.74504
## 16712 -122.4194 37.75596
## 16755 -122.3895 37.73005
## 16778 -122.4708 37.75810
## 16796 -122.4035 37.72002
## 16800 -122.3696 37.72900
## 16802 -122.4042 37.71115
## 16921 -122.3696 37.72900
## 16934 -122.4054 37.77794
## 17006 -122.4459 37.77465
## 17051 -122.4560 37.75650
## 17073 -122.4386 37.73261
## 17084 -122.3942 37.77181
## 17091 -122.4750 37.71499
## 17156 -122.4286 37.71810
## 17168 -122.4766 37.76075
## 17199 -122.3910 37.78938
## 17200 -122.4151 37.78524
## 17223 -122.4011 37.77047
## 17249 -122.4159 37.74769
## 17252 -122.4908 37.73143
## 17265 -122.4043 37.80443
## 17304 -122.4014 37.77974
## 17345 -122.4245 37.79493
## 17358 -122.4704 37.72549
## 17384 -122.3904 37.78271
## 17404 -122.4080 37.77576
## 17410 -122.3904 37.75368
## 17465 -122.5033 37.74841
## 17491 -122.4015 37.73255
## 17510 -122.4359 37.75497
## 17517 -122.4519 37.76025
## 17542 -122.4124 37.77093
## 17562 -122.4514 37.70881
## 17613 -122.3893 37.73988
## 17618 -122.4249 37.75375
## 17698 -122.4262 37.79417
## 17715 -122.4140 37.75590
## 17754 -122.4360 37.80400
## 17767 -122.4916 37.78122
## 17857 -122.4065 37.77812
## 17867 -122.4777 37.76001
## 17876 -122.5107 37.77399
## 17881 -122.4226 37.73575
## 17887 -122.4188 37.77159
## 17912 -122.3886 37.77159
## 17929 -122.4289 37.72161
## 17933 -122.4452 37.71900
## 17936 -122.4458 37.72772
## 17986 -122.4429 37.73288
## 17988 -122.4473 37.73331
## 18002 -122.4422 37.77803
## 18018 -122.4772 37.74409
## 18050 -122.4198 37.76161
## 18053 -122.4038 37.76405
## 18059 -122.4221 37.74746
## 18071 -122.4271 37.73809
## 18107 -122.4388 37.73398
## 18164 -122.4160 37.73504
## 18227 -122.4805 37.75531
## 18231 -122.4671 37.76370
## 18323 -122.4229 37.80586
## 18363 -122.3727 37.73057
## 18386 -122.4767 37.77714
## 18433 -122.4408 37.78769
## 18439 -122.4345 37.73587
## 18444 -122.4716 37.78657
## 18453 -122.3928 37.78205
## 18470 -122.4682 37.75544
## 18504 -122.5075 37.74854
## 18511 -122.3959 37.79054
## 18553 -122.3942 37.77181
## 18584 -122.4140 37.75590
## 18619 -122.4315 37.73498
## 18627 -122.4716 37.75516
## 18642 -122.4780 37.78258
## 18653 -122.4258 37.72884
## 18718 -122.4391 37.75547
## 18733 -122.4294 37.72902
## 18742 -122.4194 37.75596
## 18792 -122.4181 37.78824
## 18795 -122.3938 37.78638
## 18801 -122.4288 37.72584
## 18871 -122.4234 37.77253
## 18886 -122.4305 37.74402
## 18913 -122.4280 37.74869
## 18960 -122.4387 37.71054
## 18968 -122.4041 37.80347
## 18971 -122.4058 37.80159
## 18990 -122.4398 37.75879
## 19015 -122.4613 37.76457
## 19065 -122.4121 37.77187
## 19082 -122.4684 37.76042
## 19105 -122.4836 37.71006
## 19143 -122.4642 37.78458
## 19153 -122.4520 37.73468
## 19197 -122.4102 37.75679
## 19202 -122.4412 37.71904
## 19227 -122.3896 37.78063
## 19281 -122.4309 37.72308
## 19325 -122.4287 37.77642
## 19375 -122.4152 37.73366
## 19532 -122.4596 37.78655
## 19535 -122.5019 37.77481
## 19570 -122.4623 37.73544
## 19626 -122.4348 37.74426
## 19641 -122.3919 37.78868
## 19718 -122.4194 37.75596
## 19765 -122.4194 37.75596
## 19782 -122.3911 37.73912
## 19802 -122.4152 37.79845
## 19811 -122.3951 37.72508
## 19832 -122.4798 37.78437
## 19838 -122.4301 37.74120
## 19845 -122.4435 37.75964
## 19880 -122.4243 37.79352
## 19897 -122.3934 37.77764
## 19932 -122.4590 37.71972
## 19933 -122.4261 37.73585
## 19954 -122.4934 37.76073
## 19991 -122.4242 37.79788
## 20009 -122.4528 37.71555
## 20026 -122.3909 37.75967
## 20061 -122.4382 37.80320
## 20095 -122.3919 37.78868
## 20102 -122.4936 37.74282
## 20221 -122.4251 37.75560
## 20235 -122.4709 37.75704
## 20293 -122.4937 37.77999
## 20310 -122.4252 37.74835
## 20320 -122.4355 37.71758
## 20321 -122.4739 37.74485
## 20347 -122.5013 37.76014
## 20356 -122.4517 37.73356
## 20426 -122.4014 37.77974
## 20484 -122.4345 37.80458
## 20562 -122.3974 37.78575
## 20605 -122.4729 37.77471
## 20611 -122.3922 37.78607
## 20655 -122.4621 37.77642
## 20687 -122.4548 37.78306
## 20743 -122.4469 37.71901
## 20754 -122.4493 37.77461
## 20768 -122.3962 37.77393
## 20778 -122.4743 37.76473
## 20779 -122.4024 37.76117
## 20804 -122.5022 37.74887
## 20871 -122.3922 37.73197
## 20877 -122.4180 37.76654
## 20928 -122.4731 37.77999
## 20971 -122.4816 37.77795
## 20979 -122.4582 37.73860
## 21007 -122.3696 37.72900
## 21053 -122.4194 37.75596
## 21064 -122.4263 37.76856
## 21079 -122.4425 37.80069
## 21093 -122.4336 37.74361
## 21105 -122.4245 37.79920
## 21124 -122.4399 37.79230
## 21152 -122.4642 37.72043
## 21174 -122.4227 37.79931
## 21201 -122.4340 37.77121
## 21231 -122.4272 37.72286
## 21255 -122.4094 37.75902
## 21374 -122.4208 37.73601
## 21424 -122.3702 37.72838
## 21455 -122.4195 37.80612
## 21456 -122.3932 37.78006
## 21458 -122.4980 37.75912
## 21479 -122.4331 37.79233
## 21486 -122.3922 37.78607
## 21496 -122.4739 37.78414
## 21499 -122.4031 37.79058
## 21508 -122.4606 37.72592
## 21522 -122.4188 37.77159
## 21599 -122.4630 37.74288
## 21626 -122.4271 37.78448
## 21632 -122.4368 37.74525
## 21644 -122.4176 37.76271
## 21652 -122.3901 37.75370
## 21678 -122.3702 37.72838
## 21694 -122.4414 37.79054
## 21711 -122.4126 37.77870
## 21712 -122.3904 37.78139
## 21790 -122.4446 37.79053
## 21791 -122.5024 37.73799
## 21834 -122.4218 37.79478
## 21889 -122.4157 37.79792
## 21932 -122.4520 37.71732
## 21957 -122.4282 37.79869
## 21972 -122.4552 37.74086
## 21985 -122.4549 37.73345
## 22007 -122.4195 37.74526
## 22024 -122.5006 37.75075
## 22064 -122.4569 37.75564
## 22085 -122.4497 37.72884
## 22111 -122.3933 37.77891
## 22167 -122.5041 37.73937
## 22254 -122.3992 37.79825
## 22255 -122.4107 37.77166
## 22273 -122.4448 37.71416
## 22324 -122.3986 37.71395
## 22340 -122.4376 37.75052
## 22499 -122.4438 37.72965
## 22531 -122.4141 37.77608
## 22533 -122.4409 37.74887
## 22549 -122.4431 37.71275
## 22554 -122.4009 37.76461
## 22600 -122.4362 37.74701
## 22607 -122.4753 37.71516
## 22673 -122.4131 37.74787
## 22739 -122.3919 37.78868
## 22801 -122.4709 37.75504
## 22802 -122.4014 37.77974
## 22820 -122.4887 37.77819
## 22870 -122.4367 37.76372
## 22871 -122.4507 37.72158
## 22896 -122.4342 37.74955
## 22924 -122.4460 37.71426
## 22931 -122.4210 37.74497
## 22951 -122.4948 37.77410
## 23000 -122.4254 37.76334
## 23036 -122.4163 37.80401
## 23045 -122.3996 37.77438
## 23056 -122.4753 37.71487
## 23074 -122.3696 37.72900
## 23117 -122.4223 37.79267
## 23135 -122.3942 37.77181
## 23154 -122.4648 37.75510
## 23172 -122.4415 37.78044
## 23173 -122.4835 37.78257
## 23251 -122.4119 37.71432
## 23267 -122.4263 37.71022
## 23308 -122.4769 37.76468
## 23356 -122.4613 37.71719
## 23400 -122.4689 37.77570
## 23433 -122.4753 37.71553
## 23467 -122.3928 37.78205
## 23478 -122.4855 37.75523
## 23486 -122.4276 37.72233
## 23503 -122.4108 37.76534
## 23519 -122.4454 37.76390
## 23535 -122.4639 37.76393
## 23540 -122.4599 37.72230
## 23542 -122.4860 37.78680
## 23598 -122.3910 37.78938
## 23601 -122.3922 37.78607
## 23721 -122.4452 37.71600
## 23782 -122.4012 37.75378
## 23790 -122.5038 37.75594
## 23828 -122.3957 37.73636
## 23881 -122.3906 37.75501
## 23895 -122.4517 37.72883
## 23896 -122.4011 37.77047
## 23913 -122.4734 37.75749
## 23918 -122.4676 37.71048
## 23919 -122.4936 37.77938
## 23991 -122.4418 37.73477
## 24130 -122.4544 37.76267
## 24185 -122.4891 37.77817
## 24223 -122.4194 37.75596
## 24282 -122.4934 37.73923
## 24314 -122.4577 37.74434
## 24326 -122.3880 37.73296
## 24374 -122.3696 37.72900
## 24393 -122.4372 37.75934
## 24401 -122.4517 37.73329
## 24402 -122.4711 37.74663
## 24504 -122.4536 37.70915
## 24546 -122.4544 37.78545
## 24554 -122.4395 37.75941
## 24574 -122.3997 37.75216
## 24607 -122.4836 37.77285
## 24639 -122.4824 37.75738
## 24659 -122.3920 37.72865
## 24674 -122.4760 37.78306
## 24700 -122.4138 37.73519
## 24716 -122.4922 37.77477
## 24724 -122.4374 37.75893
## 24742 -122.4704 37.78430
## 24781 -122.4137 37.79099
## 24794 -122.4016 37.75125
## 24816 -122.4009 37.72799
## 24840 -122.4704 37.74936
## 24872 -122.4236 37.77800
## 24896 -122.4441 37.73437
## 24922 -122.4446 37.77193
## 24940 -122.4192 37.72509
## 24949 -122.4356 37.72596
## 24966 -122.3919 37.78868
## 24983 -122.4458 37.73512
## 24986 -122.4668 37.77558

str(sfhomes15) # the data frame

str(sfhomes15_sp) # the data frame

str(sfhomes15) # the data frame

## 'data.frame':    835 obs. of  19 variables:
##  $ FiscalYear      : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ SalesDate       : chr  "2015-08-21" "2015-08-13" "2015-12-29" "2015-08-26" ...
##  $ Address         : chr  "0000 2760 19TH                AV0015" "0000 0560AMISSOURI            ST0000" "0000 0718 LONG BRIDGE         ST1202" "0000 0103 SUMMIT              WY0000" ...
##  $ YearBuilt       : int  1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
##  $ NumBedrooms     : int  2 2 2 3 3 3 3 2 0 2 ...
##  $ NumBathrooms    : int  2 2 2 3 3 3 3 2 0 2 ...
##  $ NumRooms        : int  5 5 5 6 7 7 7 5 0 0 ...
##  $ NumStories      : int  0 1 1 0 1 14 2 1 0 1 ...
##  $ NumUnits        : int  1 1 1 1 1 1 1 1 0 1 ...
##  $ AreaSquareFeet  : int  1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
##  $ ImprovementValue: int  432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
##  $ LandValue       : int  432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
##  $ Neighborhood    : chr  "West of Twin Peaks" "Potrero Hill" "Mission Bay" "Lakeshore" ...
##  $ Location        : chr  "(37.7360097396496, -122.474067310226)" "(37.759197817252, -122.396516184449)" "(37.7718055392903, -122.394186974738)" "(37.7151803909173, -122.475425717555)" ...
##  $ SupeDistrict    : int  7 10 6 7 9 3 8 2 7 6 ...
##  $ totvalue        : int  865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
##  $ SalesYear       : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ lat             : num  37.7 37.8 37.8 37.7 37.8 ...
##  $ lon             : num  -122 -122 -122 -122 -122 ...

str(sfhomes15_sp) # the SPDF

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  835 obs. of  17 variables:
##   .. ..$ FiscalYear      : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   .. ..$ SalesDate       : chr [1:835] "2015-08-21" "2015-08-13" "2015-12-29" "2015-08-26" ...
##   .. ..$ Address         : chr [1:835] "0000 2760 19TH                AV0015" "0000 0560AMISSOURI            ST0000" "0000 0718 LONG BRIDGE         ST1202" "0000 0103 SUMMIT              WY0000" ...
##   .. ..$ YearBuilt       : int [1:835] 1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
##   .. ..$ NumBedrooms     : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
##   .. ..$ NumBathrooms    : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
##   .. ..$ NumRooms        : int [1:835] 5 5 5 6 7 7 7 5 0 0 ...
##   .. ..$ NumStories      : int [1:835] 0 1 1 0 1 14 2 1 0 1 ...
##   .. ..$ NumUnits        : int [1:835] 1 1 1 1 1 1 1 1 0 1 ...
##   .. ..$ AreaSquareFeet  : int [1:835] 1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
##   .. ..$ ImprovementValue: int [1:835] 432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
##   .. ..$ LandValue       : int [1:835] 432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
##   .. ..$ Neighborhood    : chr [1:835] "West of Twin Peaks" "Potrero Hill" "Mission Bay" "Lakeshore" ...
##   .. ..$ Location        : chr [1:835] "(37.7360097396496, -122.474067310226)" "(37.759197817252, -122.396516184449)" "(37.7718055392903, -122.394186974738)" "(37.7151803909173, -122.475425717555)" ...
##   .. ..$ SupeDistrict    : int [1:835] 7 10 6 7 9 3 8 2 7 6 ...
##   .. ..$ totvalue        : int [1:835] 865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
##   .. ..$ SalesYear       : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   ..@ coords.nrs : int [1:2] 19 18
##   ..@ coords     : num [1:835, 1:2] -122 -122 -122 -122 -122 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:835] "24" "35" "76" "79" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] -122.5 37.7 -122.4 37.8
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

You can see from str(sfhomes) that a SPDF object is a collection of slots or components. The key ones are:

• @data data frame of attributes that describe each location
• @coords the coordinates for each geometric object - here points
• @bbox the min and max bounding coordinates
• @proj4string the coordinate reference system defintion as a string

Explore the object in the Environment window

Review the output of each of these:

summary(sfhomes15_sp)
head(sfhomes15_sp@data)
class(sfhomes15_sp@data)

sfhomes15_sp@bbox
bbox(sfhomes15_sp)

head(sfhomes15_sp@coords)
head(sfhomes15_sp$lat)
head(sfhomes15_sp$lon)

sfhomes15_sp@proj4string
proj4string(sfhomes15_sp)

Look at sfhomes15_sp in the environment window

Are all the columns that were present in sfhomes15 also in sfhomes15_sp?

Is there a slot in sfhomes15_sp without data?

proj4string(sfhomes15_sp) # get a CRS object

## [1] NA

plot(sfhomes15_sp)  # using sp::plot

We created great maps of sfhomes point data with ggplot and ggmap.

Then we created a simple map of the SPDF with plot.

We aren't seeing the value of sp objects just yet.

rgdal is an R port of the powerful and widely used GDAL library.

It is the most commonly used R library for importing and exporting spatial data.

• OGR: for vector data: readOGR() and writeOGR()

• GDAL for raster data: readGDAL() and writeGDAL()

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.4.3

## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-5

# See what file types are supported by rgdal drivers
# ogrDrivers()$name

gdal.org

?readOGR

For more info on working with rgdal to load different types of spatial data in R see this excellent tutorial by Zev Ross.

Take a look at the file(s)

dir("data", pattern="sf_boundary")

## [1] "sf_boundary.dbf" "sf_boundary.prj" "sf_boundary.shp" "sf_boundary.shx"

sfboundary <- readOGR(dsn="data",layer="sf_boundary")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_boundary"
## with 1 features
## It has 3 fields

# or
# sfboundary <- readOGR("data","sf_boundary")
# but not
#sfboundary <- readOGR(dsn="data/",layer="sf_boundary")

What type of sp object is sfboundary? How many features are in the object?

class(sfboundary)
str(sfboundary) 
head(sfboundary@data)  

Explore the object in the Envi window

How?

plot(sfboundary)

Is it a geographic or projected CRS?

proj4string(sfboundary)

plot(sfboundary)
points(sfhomes15_sp, col="red")

Where are the points? What's wrong?

plot(sfboundary)
points(sfhomes15_sp, col="red")

Compare the CRSs, are they the same?

proj4string(sfboundary)
proj4string(sfhomes15_sp)
proj4string(sfboundary) == proj4string(sfhomes15_sp)

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] NA

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] NA

sfboundary@bbox

##         min       max
## x  542696.6  556659.9
## y 4173563.7 4185088.6

sfhomes15_sp@bbox

##            min        max
## lon -122.51072 -122.36964
## lat   37.70854   37.80612

The #1 reason…

All sp objects should have a defined CRS

If not, one needs to be assigned to the object

• This is called defining a projection

• This doesn't change the coordinates

All sp objects should have the same CRS.

• When they don't, they need to be transformed to a common CRS.

• This is called a projection transformation,
  • or projecting or reprojection.

• Projection transformation returns a new spatial object with the transformed coordinates

GIS software - or software for working with geospatial data - will have a database of definitions for thousands of Earth referenced coordinate systems

and methods for using those CRS definitions to define or transform a CRS.

So, do we need to define the CRS of either

• sfboundary or

• sfhomes15_sp?

If yes, which one? How can you tell?

We need to know

• what the appropriate CRS is
• how to define it as a CRS object
• how to assign it to an sp object

bbox(sfhomes15_sp)

##            min        max
## lon -122.51072 -122.36964
## lat   37.70854   37.80612

Use the sp function CRS() to define a coordinate reference system

• requires as input a string that defines the CRS parameters

Then, use the sp function proj4string() to assign it

proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
                               +ellps=WGS84 +datum=WGS84 +no_defs")  

The WGS84 geographic coordinate referenece system can be defined by:

1. It's paramater string

CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

1. A code that can be used to look up the parameters

CRS("+init=epsg:4326")

EPSG refers to the group that defined these

• European Petroleum Survey Group

EPSG codes are more commonly used. Can you guess why?

# use an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326") 

# or enter the parameter string
# proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
#                               +ellps=WGS84 +datum=WGS84 +no_defs")  

Proj4 is the standard library for defining and transforming map projections and coordinate reference systems.

Here is an example file of proj4 CRS definitions

Once you set the CRS you can check it with the same function.

proj4string(sfhomes15_sp)

## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

proj4string() can get or set the CRS

CRS() is used to define the CRS of an sp object

The following are equivalent

# Define the CRS with the parameter string
proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
                              +ellps=WGS84 +datum=WGS84 +no_defs")  

# Define the CRS with an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326") 

epsg codes are used as short-hand for CRSs definitions

Are they both defined? Equal?

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] FALSE

Use sp function spTransform

Requires as input:

• a sp object to transform with a defined CRS

• a target CRS

Outputs a new spatial object with coordinate data in the target CRS

• It does not change the input data

Let's transform the CRS of sfboundary to be the same as that of sfhomes15_sp.

Why not the other way around?

sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))

Do they match now?

How will we know?

proj4string(sfhomes15_sp) == proj4string(sfboundary_lonlat)

## [1] FALSE

plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")

Woo-hoo!

We can transform sp objects to the same CRS for mapping and spatial analysis (we'll get there!)

Use writeOGR to save sfboundary_lonlat to a new shapefile

See ?writeOGR for help

# write transformed data to a new shapefile 
writeOGR(sfboundary_lonlat, 
          dsn = "data", 
          layer = "sfbounary_lonlat", 
          driver="ESRI Shapefile")

# is it there?
dir("data")

We want all data in the same CRS

Which one is best?

Geographic CRSs

• 4326 Geographic, WGS84 (default for lon/lat)

• 4269 Geographic, NAD83 (USA Fed agencies like Census)

Projected CRSs

• 5070 USA Contiguous Albers Equal Area Conic

• 3310 CA ALbers Equal Area

• 26910 UTM Zone 10, NAD83 (Northern Cal)

• 3857 Web Mercator (web maps)

See http://spatialreference.org/

Use this site to find EPSG codes and proj4 CRS strings

1. Use coordinates to convert the landmarks data frame to a SpatialPointsDataFrame
2. Use proj4string to define it's CRS as 3857 Web Mercator
3. Use spTransform to transform the landmarks data to geographic coords (WGS84)
4. Add the landmarks to the last map we created.

# Convert the DF to a SPDF
coordinates(landmarks) <-c("X","Y")
# Define the CRS with an EPSG code for Web mercator
proj4string(landmarks) <- CRS("+init=epsg:3857") 
# Transform the CRS to WGS84
landmarks_lonlat <- spTransform(landmarks, CRS("+init=epsg:4326"))
# map it
plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")

1. Use readOGR to load the shapefile "data/sf_highways.shp"
2. Check its data type with class()
3. Check its CRS with proj4string()
4. Use spTransform if need to project it to WGS84
5. Recreate the last map and add the highways as black lines

highways <- readOGR(dsn="data", layer="sf_highways")
class(highways)
proj4string(highways)
highways_lonlat <- spTransform(highways, CRS("+init=epsg:4326"))

plot(sfboundary_lonlat)
lines(highways_lonlat, col="black")
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_highways"
## with 246 features
## It has 5 fields

## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

# QUESTIONS?

base::plot, ggplot, ggmap for geospatial data in data frames

• great for creating maps given these types of data

sp::plot for sp objects

• meh, but great for a quick look at spatial data

Use it to create quick maps

Can be used to create great maps

• BUT complex, non-intuitive syntax = long ugly code

See examples in sp Gallery: Plotting maps with sp

# map of the sfhomes data by totalvaue
spplot(sfhomes15_sp,"totvalue")

tmap stands for thematic map

Great maps with less code than the alternatives

Syntax should be familar to ggplot2 users, but simpler

Relatively easy to create interactive maps

• ?tmap
• vignette("tmap-nutshell")
• vignette("tmap-modes")
• tmap in a Nutshell

Load the library

library(tmap)

## Warning: package 'tmap' was built under R version 3.4.3

qtm(sfhomes15_sp)

tmap_mode("view")

## tmap mode set to interactive viewing

qtm(sfhomes15_sp)

tmap_mode("plot")

## tmap mode set to plotting

ttm()       # Toggle Tmap mode between "interactive" and "plot"

last_map()  # display the last map 

tmap_mode("plot")

## tmap mode set to plotting

tmap's flexibility comes in how it intitively allows you to layer spatial data and style the layers by data attributes

Use tm_shape(<sp_object>) to specifiy a geospatial data layer

Add + tm_<element>(...) to style the layer by data values

…and other options for creating a publication ready map

?tmap_shape

?tmap_element

• ?tm_polygons (tm_fill, tm_borders)

• ?tm_symbols (tm_dots, etc…)

• tmap in a Nutshell

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black")

tm_shape(sfhomes15_sp) + 
  tm_dots(col="red", size=.25)

tm_shape(highways_lonlat) + 
  tm_lines(col="black")

## Mapping values

tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25)    # column names must be quoted

tm_dots are a type of tm_symbols()

See ?tm_symbols for other types of point symbols.

tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25)  # columnn names **must** be quoted!

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25) 

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25, title = "San Francisco Property Values (2015)") + 
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

Add the landmarks to the previous tmap

Feel free to change the colors and title etc.

What does the landmarks layer tell you about the way tmap handles CRSs?

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25, 
          title = "San Francisco Property Values (2015)") + 
tm_shape(landmarks) +
  tm_markers(col="black", size=.4, text="name",text.just = "left") +
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

Make the previous tmap interactive

ttm()       # Toggle Tmap mode between "interactive" and "plot"
last_map()  # display the last map 

## tmap mode set to interactive viewing

popup.vars = c("SalesYear","totvalue","NumBedrooms", "NumBathrooms","AreaSquareFeet")

tm_shape(sfboundary_lonlat) + 
  tm_polygons(col="beige", border.col="black") +
tm_shape(highways_lonlat) + 
  tm_lines(col="black") +
tm_shape(sfhomes15_sp) + 
  tm_dots(col="totvalue", size=.25, 
          title = "San Francisco Property Values (2015)",
          popup.vars=c("SalesYear","totvalue","NumBedrooms",
                       "NumBathrooms","AreaSquareFeet")) + 
tm_shape(landmarks) +
  tm_markers(col="black", size=.4, text="name",text.just = "left") +
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

Save it as an R object

map1 <- last_map()
map1 # then display it

Try this and then take a look at your files

save_tmap(map1, "sf_properties.png", height=6) # Static image file
save_tmap(map1, "sf_properties.html") # interactive web map

## Map saved to /Users/pattyf/Documents/Dlab/workshops/2018/rgeo_may2018/r-geospatial-workshop/sf_properties.png

## Resolution: 2318.548 by 1800 pixels

## Size: 7.728492 by 6 inches (300 dpi)

## Interactive map saved to /Users/pattyf/Documents/Dlab/workshops/2018/rgeo_may2018/r-geospatial-workshop/sf_properties.html

Many ways to do this.

You can share you map online by publishing it to RPubs.

• You need to have an RPubs account to make that work.

1. Enter the name of your tmap object (map1) or your tmap code in the console

2. In the Viewer window, click on the Publish icon.

tmap_mode("interactive")
map1 #  in console

http://rpubs.com/pfrontiera/358110

# write transformed data objects to a new shapefile 
coordinates(bart) <- c("X","Y")
proj4string(bart) <- CRS(proj4string(sfboundary_lonlat))

save_layers <- c("sfhomes15_sp", "sfboundary_lonlat", "highways_lonlat", "bart", "landmarks")
for (i in save_layers){
  print(paste("Saving [",i,"] to file..."))
  writeOGR(get(i), 
          dsn = "saved_data", 
          layer = i, 
          driver="ESRI Shapefile")
}
# is it there?
dir("saved_data")

• Geospatial Data in R
• CSV > Data Frame > ggplot/ggmap
• sp library - spatial objects and methods in R
• GDAL libary - readOGR, writeOGR
• Convert Data Frame to SPDF
• CRS definitions and transformations
• Map Overlays
• Geospatial Maps with tmap